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Abstract 

The first- and second-order correlation functions of trapped, interacting Bose-Einstein conden- 
sates are investigated numerically on a many-body level from first principles. Correlations in real 
space and momentum space are treated. The coherence properties are analyzed. The results are 



obtained by solving the many-body Schrodinger equation. It is shown in an example how many- 

(N 

body effects can be induced by the trap geometry. A generic fragmentation scenario of a condensate 
is considered. The correlation functions are discussed along a pathway from a single condensate to 

-5 

a fragmented condensate. It is shown that strong correlations can arise from the geometry of the 
trap, even at weak interaction strengths. The natural orbitals and natural geminals of the system 



are obtained and discussed. It is shown how the fragmentation of the condensate can be under- 
stood in terms of its natural geminals. The many-body results are compared to those of mean-field 

O 

theory. The best solution within mean-field theory is obtained. The limits in which mean-field the- 

i— I ories are valid are determined. In these limits the behavior of the correlation functions is explained 
> 

I" — within an analytical model. 
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I. INTRODUCTION 



The computation of correlation functions in interacting quantum many-body systems is 
a challenging problem of contemporary physics. Correlations between particles can exist 
in time, in real space or in momentum space. Of course, all combinations of the above 
three cases are possible. Since the first experimental realization of Bose-Einstein conden- 
sates (BECs) in ultracold atomic gases [1-3], great experimental and theoretical progress 
has been made in the determination of the coherence and the correlation functions of Bose- 
condensed systems. Over the years experiments have measured more and more accurately 
first, second and to some extent even third order correlations of trapped BECs, see [4-11]. 
Theoretically, the correlation functions of trapped, interacting BECs have been investigated 
in numerous works, see e.g. [12-21]. While analytical approaches from first principles are 
usually restricted to treat homogeneous gases without any trapping potential, numerical 
methods can overcome this restriction. It is important to note that the shape of the trap- 
ping potential can have a substantial impact on the properties of the many-body system. 
This is particularly true for issues concerning condensation [22] and fragmentation of Bose 
systems [23]. For example, the ground state of weakly interacting condensates in harmonic 
traps is almost fully condensed, while the ground state of double-well potentials can be 
fragmented or condensed, depending on the height of the barrier, the number of particles 
and the interaction strength [24-29] . In this work we investigate first- and second-order cor- 
relations of trapped, interacting condensates and their coherence properties depending on 
the trap geometry from first principles. Our results are obtained by solving the many-body 
Schrodinger equation of the interacting system numerically. From this many-body solution 
we extract the first- and second-order reduced density matrices which allow us to compute 
all real and momentum space first- and second-order correlations and in particular the frag- 
mentation of the condensate. For illustration purposes we consider a stationary system in 
the ground state to show how many-body effects can become dominant when the trap geom- 
etry is varied. As a numerical method to solve the interacting many-body problem we use 
the recently developed multiconfigurational time- dependent Hartree for bosons (MCTDHB), 
which propagates a given many-body state in time [30-32]. By propagation in imaginary 
time it allows us to investigate the ground state and other stationary states. Alternatively, 
one can use the stationary multiconfigurational Hartree for bosons to compute these states 
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[29]. In order to identify true many-body effects in the correlation functions, we compare our 
many-body results with those based on mean-field approaches. More specifically, we com- 
pute the energetically-lowest mean-field solution of the same system for comparison. This is 
the best approximation to the true many-body wave-function within mean-field theory and, 
thereby, allows us to pinpoint the limits of mean-field theory. A general method to compute 
this best mean-field solution has been developed in our group [25-28, 33]. For completeness 
we compare the many-body results also to the widely used Gross-Pitaevskii mean-field so- 
lution. In order to understand first- and second-order correlations in an intuitive way, we 
develop an analytical mean-field model which explains the general structure of our results 
in those regions where many-body effects can be safely neglected. 

This paper is organized as follows. In Sec. II we review basic facts about reduced density 
matrices, correlation functions and coherence. In Sec. Ill we give a brief introduction to 
the numerical many-body and mean-field methods that we use. In Sec. IV we introduce a 
generic one-dimensional model system that we solve. In particular, we identify mean-field 
and many-body regimes of this system. We explain the transition from a condensed to a 
fragmented state in terms of the natural geminals of the system. In Sees. V and VI we 
present our many-body results for the first- and second-order correlation functions and the 
coherence properties of the model system. We compare our many-body results with those 
obtained by using the best mean-field and the Gross-Pitaevskii approximation. The structure 
of the correlation functions and the coherence are explained within an analytical mean-field 
model in the limits where mean-field theory is applicable. 

II. BASIC DEFINITIONS 

We consider a given wave function \l/(ri, . . . ,r^;t) of N identical, spinless bosons with 
spatial coordinates in D dimensions. The p-th order reduced density matrix (RDM), is 
defined by 

JV' f 

p w (ri, . . . , r p \r[, . . . , r' p ;t) = _ j tf(ri, . . . , r p , r p+1 , . . . , r N ; t) 

x^*(r'i, ■ • • , rj,, r p+1 , . . . , r N ; t)dr p+1 . . . dr N , (1) 
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where the wave function is assumed to be normalized {^(t)\^(t)) = 1. Equivalent to Eq. (1), 
p( p \ri, . . . , r p \r[, . . . , r^; t) can be regarded as the kernel of the operator 

p (p) = jj/^^N-AMt))m)\] (2) 

in Hilbert space, where Trjv_ p [-] specifies taking the partial trace over N — p particles. Since 
the wave function is symmetric in its coordinates, it does not matter over which particles 
the trace is taken. In what follows, we add |^) as an additional subscript if a result is only 
valid for states |^) of a particular form. 

The diagonal p (j '\r 1 , . . . , r p |r i, . . . , r p ; t) is the p-particle probability distribution at time 
t multiplied by N\/(N— p)\. The p-th order RDM p^ can be expanded in its eigenfunctions, 
leading to the representation 

p^( ri , . . . , r p K, . . . , r; ; t) = Yl ^K (p) (n, ...,r p , Vi, • • • , t). (3) 

i 

Here, nf\t) denotes the i-th eigenvalue of the p-th order RDM and otf\vi,. . . ,r p ,t) the 
corresponding eigenfunction. The eigenfunctions are known as natural p-functions and the 
eigenvalues as natural occupations. For p = 1 and p = 2 the eigenfunctions are also known 
as natural orbitals and natural geminals, respectively. We order the eigenvalues n[ p \t) for 
every p non-increasingly, such that nf\t) denotes the largest eigenvalue of the p-th order 
RDM. The normalization of the many-body wave function and Eqs. (1) and (3) put the 
restriction 

on the eigenvalues of the p-th order RDM. Thus the largest eigenvalue nf\t) is bounded 
from above by [34, 35] 

Lower bounds on nf\t) can be derived, relating RDMs of different order [36, 37]. In 
particular, for the case p = 2 it can be shown that [36] 

nf\t)>n { l\t)[n^\t)-l}. (6) 

It is a well known fact that the natural orbitals of a symmetric (or antisymmetric) function 
constitute a sufficient one-particle basis to expand \1/ and the eigenfunctions of the RDMs 



for all p [34]. It is therefore possible to construct the natural geminals in the basis of 
one-particle functions spanned by the natural orbitals. 

The determination of accurate bounds on eigenvalues of RDMs is an active field of research 
[34, 35] since it is possible to express the exact energy expectation value of a quantum system 
of identical particles interacting via two-body interactions by an expression involving only 

(2) (2) 

the natural geminals, a\ (ri, r 2 , t), and their occupations, n\ (t). For a general Hamiltonian 

N N 

If = + £ W ( r '- r i)> (7) 

i=l i<j 

consisting of one-body operators hfa) and two-body operators Wfa — r^), the expectation 
value of the energy E can be expressed by making use of the time-dependent natural geminals 
af , \xi,x 2 ,t) and following Refs. [34, 35] through the equation: 

'h(n) + h(r 2 ) 



N- 1 



+ W(r l -r 2 ) 



afVi,^)- (8) 



Note that the many-body wave function does not appear explicitly in Eq. (8). We will not go 
any further into the details of these approaches to many-body physics and refer the reader 
to the literature [34, 35]. 

The natural orbitals also serve to define Bose-Einstein condensation in interacting sys- 
tems. According to Penrose and Onsager [22], a system of identical bosons is said to be 
condensed, if the largest eigenvalue of the first-order RDM is of the order of the number of 
particles in the system, = O(N). If more than one eigenvalue of the first-order RDM is 
of the order of the number of particles, the condensate is said to be fragmented, according 
to Nozieres and Saint James [23]. 

Equivalent to Eq. (1), the p-th order RDM can be expressed through field operators as 

p«(n, . . . , r p |ri, . . . , r' p , t) = ^(t)\^(r[) . . . * t (r;)*(r p ) . . . *(n)|*(t)>, (9) 

where the Schrodinger field operators satisfy the usual bosonic commutation relations 

^(r),*^')] =S(r-r'), [*(r), *(r')] = 0. (10) 

The representation given in Eq. (9) shows that the p-th order RDM is identical to the p-th 
order correlation function at equal times [12, 38]. 
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In order to discuss correlations not only in real space, but also in momentum space, we 
define the Fourier transform of a function /(r 1? . . . , r p ) of p D-dimensional coordinates by 

f(k 1 ,...,k p ) = ^-^ 2 |^re-^r- k ' r '/(ri,...,r p ). (11) 

By applying the Fourier transform, Eq. (11), to the coordinates ri, . . . , r p and r[, . . . , r' p of 
the natural p-functions a^\ri, . . . , r p , t) and af^i^'i, ■ ■ ■ , r ' p ,t) in Eq. (3), one arrives at the 
momentum space representation of pW: 

pW(k!, . . . ,k p |k;, . . . ,k' p ,t) = ^^(t)^^, . . . ,k p ,tK (p) Vi, • • • ,k;,t). (12) 

i 

The diagonal p^(ki, . . . , k p |ki, . . . , k p ; t) in momentum space is the p-particle momen- 
tum distribution at time t, multiplied by N\/(N — p)\. It can be shown that the 
p-particle momentum distribution at large momenta is dominated by contributions of 
p( p \ri, . . . ,r p \r[, . . . ,r' p ;t) close to the diagonal, i.e. r- for % = I,..., p. Sim- 

ilarly, the p-particle distribution at low momenta is dominated by the behavior of 
p^(ri, . . . , r p \r' 1: . . . , r' p ; t) on the off-diagonal at large distances between and r-. See 
Appendix A for more details. 

Apart from the p-particle distributions themselves, either in real space or in momentum 
space, it is also of great interest to compare the p-particle probabilities to their respective 
one-particle probabilities. Thereby, it becomes possible to identify effects that are due to 
the quantum statistics of the particles. The normalized p-th order correlation function at 
time t is defined by [38] 

(p)/ / / ,\ p^(ri, . . . , rpjr^, . . . ,r' p ;t) 

9 (r '-" r ""-- r ' ;,) " VE^*P»P) ( ) 

and is the key quantity in the definition of spatial coherence. Full spatial p-th order coherence 
is obtained if p( n \r 1: . . . r p \r' 1: . . . , r' p ; t) factorizes for all n < p into a product of one complex 
valued function S(r,t) of the form 

. . . , r p |r;, . . . , r' p , t) = £*(r[, t) ■ ■ ■ £*(v' n , t)£(v n , t) ■ ■ ■ S(r 1} t). (14) 

In this case 

|^(r;,...,r;, ri ,...,r p ;i)| = l (15) 

for all n < p. Otherwise, the state \^{t)) is only partially coherent. Full coherence in a 
system with a definite number of particles iV can only be obtained for p — 1 [39]. However, 

6 



when the particle number N is large, p-th order coherence can be obtained up to corrections 
D(l/N), at least for p < N [39]. 

The diagonal of the normalized p-th order correlation function gW(ri, . . . , r p , r 1; . . . , r p ; t) 
gives a measure for the degree of p-th order coherence. For values 

g(p\ri, . . . ,r p ,Ti, . . . ,r p ;t) > 1 (< 1) the detection probabilities at positions r±, . . . ,r p are 
correlated (anticorrelated). 

Note that if Eq. (14) holds in real space, it must also hold in momentum space, as can 
be seen by Fourier transforming each of the 2n variables in Eq. (14). It is therefore possible 
to define the normalized p-th order correlation function in momentum space by 

S (k„...,k p ,k 1 ,...,k 5 , t ) - v / nLip(1)(ki | ki;tyi , (k ;| k:; ^ ( 16 » 

The diagonal of Eq. (16), g^(ki, . . . , k p , ki, . . . , k p ; t), expresses the tendency of p momenta 
to be measured simultaneously. For values g^(ki, . . . , k p , ki, . . . , k p ; t) > 1 (< 1) the de- 
tection probabilities of momenta ki, . . . , k p are correlated (anticorrelated). The p-th order 
momentum distribution p^(k 1; . . . , k p |k 1; . . . , k p ; t) depends on the entire p-th order RDM 
p( p )(ri, . . . , Tplr^, . . . ,r' ;t), see Appendix A. Thus, ^^(ki, . . . , k p , k 1; . . . , k p ; t) provides in- 
formation about the coherence of \^(t)) which is not contained in g( p \ri, . . . , r p , r 1 , . . . , r p ; t). 

In Young double slit experiments using non-interacting bosons ^^(r^, ri; to) I = 1 ensures 
the maximal fringe visibility of an interference pattern. Here, to is the time of release from 
the slits. See, e.g. [40], for an experiment using Bose-Einstein condensates. However, 
if interactions during the expansion behind the slit are not negligible, there is no simple 
relation between the fringe visibility and the wave function ^(i - !, . . . ,r^,t ) at the time of 
release from the slits. In other words, the interaction between the particles can modify the 
observed interference pattern [41-43]. 

In order to determine the degree of coherence of a given system, it is necessary to quan- 
tify how well Eq. (14) is satisfied. A visualization of the degree of coherence is highly 
desirable, as it helps to understand the coherence limiting factors in an intuitive man- 
ner. Already for one-dimensional systems, D — 1, the normalized first-order correlation 
function at time t, g^(r' 1: ri,t), is a complex function of two variables and cannot be 
visualized in a single plot. It is, therefore, necessary to consider quantities that sample 
parts of g^(r' 1: . . . , r p , ri, . . . , r p ; t). In one-dimensional systems j^^^r^, ri; t)\ 2 can be rep- 
resented as a two-dimensional plot and gives a measure for the degree of first-order coher- 



ence. Similarly, g^(r 1: r 2 , r 1? r 2 ; t) and (^(k]., k 2 , k 1; k 2 ; t) are real and can be represented 
as two-dimensional plots, if D — 1. In Sees. V and VI we will visualize the degree of first- 
and second-order coherence of a one-dimensional system, defined in Sec. IV, by means of 
(^(r'l, ri; t)\ 2 , 3 (2) (ri, r 2 , n, r 2 ; i) and ff ( 2 )(ki, k 2 , k 1? k 2 ; t). 

III. NUMERICAL METHODS 

The main goal of this work is to investigate exactly the behavior of first- and second-order 
correlation functions in interacting many-body systems. This requires the computation of 
the exact many-body wave function which is generally a difficult problem to solve. In some 
cases, when the general form of the wave function is known a priori, an exact solution can be 
obtained, either by solving transcendental equations or by exploiting mapping theorems, see 
e.g. [18, 44-52]. However, in general it is necessary to solve the full many-body Schrodinger 
equation numerically in an efficient way. In Sec. Ill A we give a brief account of the numerical 
method MCTDHB to solve the interacting many-boson problem. In order to find out to 
which extent mean-field methods are applicable to bosonic, interacting many-body systems, 
we compare our many-body results with those based on mean-field approaches, namely the 
commonly used Gross-Pitaevskii mean-field [53-55] and the best mean-field (BMF), which we 
describe briefly in Sec. Ill B. It is beyond the scope of this work to explain either MCTDHB 
or BMF in detail and we refer the reader to Refs. [30-32] and Refs. [25-28, 33] for more 
detailed explanations of MCTDHB and BMF, respectively. 

A. The many-body wave function 

The exact wave function of an interacting iV-boson problem can always be expanded 
in any complete set of permanents of N particles. Each of the permanents is constructed 
from a complete set of single-particle functions which are commonly referred to as orbitals. 
Practical computations can never be carried out in complete basis sets and therefore, it is 
crucial to cut the basis set carefully. 

Our starting point is the Schrodinger picture field operator *(r) satisfying the usual 
bosonic commutation relations Eqs. (10). It is convenient to expand the field operator in a 
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complete set of time- dependent orthonormal orbitals, 

*(r) = ]T4(t)0 fc (r,t), (17) 

k 

where the time-dependent annihilation and creation operators obey the usual commutation 
relations dk(t)c^(t) — c^(t)ck(t) = Skj for bosons at any time. Note that it is not necessary 
to specify the shape of the orbitals at this point. 

The many-body Hamiltonian (7) is standardly written in second quantized form as 

H = 4^*9 + \ Yl clclciCgWksqi, (18) 

k,q k,s,l,q 

where the matrix elements of the one-body Hamiltonian h(r) and two-body interaction 
potential W(r — r') are given by 



hkqit) = j (f>* k (r,t)h(r)(f) q (r,t)dr, 
W ksql (t) = yy^(r,O0:(r , ,*Wr-O^(r,^Kr , ,t)drdr'. (19) 

The ansatz for the many-body wave function in MCTDHB is taken as a linear 

combination of time-dependent permanents 

= Y C n{t) |ni, n 2 , . . . , n M ; t) , 

n 

1 / j. \ ni / , \ ni / , \ n M 

\n 1 ,n 2 ,...,n M ;t) = -j ===^ = ^ {c[(t)J [cl(t)J ■■■yc J M (t)J \vac) , (20) 

where \ni, ■ ■ ■ , um', t) is assembled from the time-dependent orbitals above. The sum- 
mation in (20) runs over all ( N+ ^~ 1 ) permanents generated by distributing iV bosons 
over M orbitals. We collect the occupations in the vector n = (ni, n 2 , ■ ■ ■ , um), where 
n\ + n 2 + . . . + n M = N. Of course, if M goes to infinity then the ansatz (20) for the wave 
function becomes exact since the set of permanents \ni,n 2 , . . . ,riM',t) spans the complete 
iV-particle Hilbert space. In practical computations we have to restrict the number M of 
orbitals from which the permanents \n\,n 2 , . . . ,riM',t) are assembled. By substituting the 
many-body ansatz (20) into the action functional of the time-dependent Schrodinger equa- 
tion, it is possible to derive a coupled set of equations of motion containing the coefficients 
Ca(t) and the set of time-dependent orbitals (j> k (r,t). The equations are obtained by re- 
quiring the stationarity of the action functional with respect to variations of the coefficients 



Cft(t) and the set of time-dependent orbitals (f>k(r,t). These coupled equations have to be 
solved simultaneously, leading to an efficient wave package propagation method for bosons 
[30-32] . At first sight it might seem to be an unnecessary complication to allow the orbitals 
0fc(r, t) to depend on time. However, this additional degree of freedom allows both, the basis 
of one-particle functions </>fc(r, t) and the coefficients Ca(t) to be variationally optimal at any 
time. Note that this is fundamentally different to a multi-mode ansatz with fixed orbitals 
in which the quality of the chosen basis set may deteriorate as the system evolves in time. 

In order to investigate stationary properties of Bose-Einstein condensates, we use a many- 
body relaxation method. By propagating a given initial guess in imaginary time with MCT- 
DHB, the system relaxes to the ground state, which allows us to treat stationary systems 
as well. A necessary requirement for this procedure to work is that the initial guess has 
non-zero overlap with the ground state. The variational principle ensures that the set of 
orbitals (pk is variationally optimal, in the sense that the lowest ground state energy within 
the Hilbert subspace of N bosons distributed over M orbitals is obtained, see in this respect 
Ref. [29]. We will not go any further into the details of MCTDHB and refer the reader to 
the literature [30-32] . 

B. Best mean- field 

The exact many-body wave function of a bosonic system of N particles can always be 
expanded in an infinite weighted sum over any complete set of permanents of iV particles. 
In mean-field theory the exact many-body wave function is approximated by a single per- 
manent. This single permanent is built from a number M < N of orthogonal orbitals in 
which the N bosons reside. In the field of Bose-Einstein condensates one particular mean- 
field, the Gross-Pitaevskii (GP) mean-field, has proven to be very successful. In analogy 
to non-interacting BECs, in GP theory it is assumed that the many-body wave function 
is given by a single permanent in which all particles reside in one orbital, i.e. M — 1. A 
minimization of the energy functional with the GP ansatz wave function leads to the famous 
Gross-Pitaevskii equation [53-55] . The solution of the GP equation yields the single orbital 
from which the GP mean-field permanent is constructed. 

However, it has been shown [25-28, 33] that the GP mean-field is not always the 
energetically-lowest mean-field solution. The assumption that all particles occupy the same 
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orbital is too restrictive. Especially in multi-well trapping geometries the energetically-lowest 
mean-field solution can be fragmented [25-28, 33], see also Sec. II. 

In order to obtain the energetically-lowest mean-field solution, it is necessary that the 
ansatz for the wave function is of the most general mean-field form. Due to the variational 
principle, the minimization of the respective energy functional with respect to all parameters 
of the ansatz wave function will then give the best solution within mean-field theory. It is 
therefore legitimate to call this mean-field solution the best mean field (BMF). A procedure 
to obtain the best mean-field (BMF) solution numerically has been developed recently [25- 
28, 33]. 

In the best mean- field approach the ansatz for the wave function |\&) is taken as a single 
permanent of N bosons distributed over M time-independent orthonormal orbitals 0fc(r): 

|^) = |m,n 2 ,...,n M ). (21) 

Using this ansatz for the wave function, the energy functional is minimized by a variation 
over the number of orbitals M, the occupation numbers rii and the orbitals 0fe(r) themselves 
[25, 28]. The variation leads to a set of coupled non-linear equations that have to be solved to 
obtain the BMF solution. Thereby, the energetically most favourable permanent is selected 
to approximate the true many-body wave function. The GP mean-field is contained in the 
BMF ansatz as can be seen by restricting the number of orbitals to M — 1. 

IV. A MODEL AND ITS PHYSICS 

In order to examine correlation functions of Bose-condensed systems, we now turn to 
a specific example. For simplicity we work in one dimension, D — 1, and henceforth we 
substitute r = x and k = k. We will study the correlation functions of N — 1000 repulsively 
interacting bosons in a double-well trap at various barrier heights. The dynamics of a similar 
system has been investigated recently in the context of a dynamically raised barrier [30]. 
In order to isolate physical effects that are due to the trapping geometry and not due to 
dynamical parameters such as the rate at which the barrier is raised, etc., we restrict our 
discussion to the ground state at different barrier heights. The restriction to a stationary 
state allows us to omit the time argument in all physical quantities from now on. Double- well 
systems have the interesting property that depending on the height of the barrier and/or the 
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interaction strength, the ground state undergoes a transition from a single to a fragmented 
condensate [24, 26, 28, 29]. We shall show how this transition from a condensed state to a 
fragmented condensate manifests itself in the correlation functions. 
We work with a dimensionless Hamiltonian of the form 

N ( l d 2 \ N 

H = Y.{-2d^ + V ^ Xi) ) +*o£; (22) 

i=l ^ * ' i<j 

to solve the stationary Schrodinger equation H^f = E^> . All quantities in Eq. (22) are 
dimensionless and the connection to a dimensional Hamiltonian 

N / j-2 a2 \ N 

H = E ( " + ^) + X ° E - **) (23) 

i=l \ i / j<j 

is made by the relations Xi = £i/L, where L is a length scale, V(xi) = ^j^V^Lxi), A = 
^Ao, • • • , xn) = &(x±L, . . . ,xnL) and E = E 11 ^-. As an external potential we choose 

a harmonic trap with an additional central barrier of variable height V(xi) = \x\ + Ae^^\ 
where A is the height of the potential barrier and a = 2 a fixed width. For the strength of 
the dimensionless interparticle interaction we choose A = 0.01. In the computations using 
MCTDHB we restrict the number of orbitals to M = 2, yielding a total of = 1001 

permanents. 



A. Condensed state 



We begin with a discussion of the ground state energy as a function of the barrier height. 
The ground state energy per particle of the many-body solution, E MC tdhb/N (blue), is 
shown in Fig. 1 (top). Emctdhb/N increases with the height of the central barrier. The 
energy differences per particle of the many-body and the BMF solution with respect to the 
GP solution, (Emctdhb — E G p)/N (blue) and (E B mf — E GP )/N (red), are shown in the 
inset of Fig. 1 (top). The energy difference (Emctdhb — E GP )/N is negative, because the 
interacting system can lower its energy by depleting the condensate. At low barrier heights 
the GP mean-field is the best mean-field and thus (Ebmf ~ E G p)/N = 0. A comparison of 
the energy scales of Fig. 1 and its inset reveals that the energy of the many-body solution, 
the BMF solution and the GP solution are very close at all barrier heights. The nature of 
the many-body ground state at different barrier heights varies nevertheless very strongly, as 
we shall show below. 
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Fig. 1 (middle) shows the occupations n^/N and n^/N of the first and second natural 
orbitals of the many-body wave function as a function of the barrier height, computed with 
MCTDHB. The largest eigenvalue of the first-order RDM, is only restricted by Eq. (5) 
and can therefore take on any value between and N. The dashed lines indicate these upper 
and lower bounds on n^{\ At low barrier heights only one natural orbital is significantly 
occupied. Therefore, we refer to the parameter range < A < 13 as the condensed regime, 
in accordance with the definition of Penrose and Onsager [22]. The occupation of the second 
natural orbital is due to the two-body interaction between the particles. However, it remains 
below 1% for all values of the barrier height A < 13 and is even below 1% at A = 0. 

Since n{ ps N in the condensed regime, the upper and the lower bounds, Eqs. (5) and 

(2) 

(6), on the largest eigenvalue of the second-order RDM, n\ , are almost identical. Therefore, 
nf* is constrained to take on a value very close to N(N — 1). Consequently, there can be only 
one significantly occupied natural geminal. This is confirmed in Fig. 1 (bottom), where the 
natural geminal occupations are shown as a function of the barrier height. For the purpose 
of describing first- and second-order correlations it is therefore legitimate to approximate 
the many-body wave function in this regime by a single permanent |iV, 0) in which all N 
bosons occupy the first natural orbital a^\xi). 

The first column of Fig. 2 shows the first (red) and the second (blue) natural orbitals 
of the many-body solution at barrier heights A = 0, 13, 19, 24, from top to bottom. The 
first and the second natural orbitals are symmetric and antisymmetric about the origin, 
respectively. At A = the first natural orbital, a^fai), takes on the shape of a broadened 
Gaussian, reflecting the repulsive interaction between the particles. The second natural 
orbital, a^\x±), has a higher kinetic energy than the first one due to the node at the center 
of the trap. Additionally, the second natural orbital forces the particles to occupy regions 
of the trap where the trapping potential is higher. There is an energy gap between the 
one-particle energies of the first and second natural orbital. The occupation of the second 
natural orbital is therefore very small in the purely harmonic trap at the chosen interaction 
strength. 

As the barrier height is varied from A = up to A = 13, the natural orbitals deform to fit 
the new shape of the external potential. The central peak of the first natural orbital splits 
into two maxima which become localized at positions X\ = ±d/2, where d is the distance 
between the wells of the external potential. 
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At the center of the trap, where the barrier is raised, the first natural orbital develops 
a local minimum in order to minimize the potential energy. The second natural orbital on 
the other hand has a node at the center of the trap at any barrier height. Its maximum and 
minimum are localized at the minima of the external potential. As the barrier is raised, the 
energy gap between the first two natural orbitals decreases. However, the increase of the 
depletion of the condensate from /N < l%o at A = to n^/N pa 1% at A — 13 cannot 
be explained in this single particle picture. On a single particle level the ground state would 
be fully condensed at any finite barrier height, i.e. = 0. The reason for the observed 
increase in the depletion lies in the fact that for repulsively interacting many-boson systems 
in multi-well setups it becomes energetically more favourable to fragment as the barrier 
between the wells is raised [24, 26, 28, 29, 33]. The increase in energy which results from the 
occupation of orbitals with a higher one-particle energy can be outweighed by a decrease 
in interaction energy. This effect becomes dominant at barrier heights above A = 13, see 
Sees. IV B and IV C. 

The second to fourth column of Fig. 2 show (from left to right) the first three natural 
geminals a\ (xi,X2) at the same barrier heights as above. By expanding the natural gem- 
inals in the one-particle basis of natural orbitals one finds that the natural geminals in the 
condensed regime are approximately given by symmetrized products of the natural orbitals: 

|«S 2) ) = |2,0), |4 2) ) = |1,1), |4 2) ) = |0,2), (24) 

where |mi,m 2 ) denotes a state with mi particles in the first and m 2 particles in the second 
natural orbital. Only the first natural geminal, a\ (xi,X2), is significantly occupied in the 
condensed regime. Due to the two-body interaction between the particles there is a small 
occupation of the second and third natural geminal. However, at low barrier heights their 
occupation is largely suppressed, due to the gap between the single particle energies of the 
first and second natural orbital. Since the geminals a 2 2 \xi,X2) and ot£\x\,X2) contain the 
second natural orbital in their expansion, see Eq. (24), their occupation increases the total 
energy at low barrier heights. 

In the equation for the energy expectation value, Eq. (8), the only substantially contribut- 
ing natural geminal is af\xi, x 2 ). The shape of a^\xi, x 2 ) is particularly interesting. It has 
four maxima of similar height, located at positions x\ = Xi = ±d/2 and x\ = — x 2 = ±d/2, 
see the second panel in the second row of Fig. 2. Since af \xi,x 2 ) has peaks on the diagonal 
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at x\ = X2 = ±d/2, it contributes to both, the one-particle part and the interaction part of 
the energy. 

In contrast to the first natural geminal, a 2 (xi,x 2 ) and a% (xi,x 2 ) both exhibit node 
lines going through the region where the central barrier is raised. As the energy gap between 
the single particle energies of the natural orbitals a{\xi,X2) and a^\xi,X2) decreases, so 
does the energy gap between the natural geminals. Similar to the discussion of the natural 
orbital occupations above, this argument in terms of an energy gap does not explain the 
increase of the occupation of the second and third natural geminal when the barrier is raised. 
Without interactions the occupation numbers of all but the first natural geminal would be 
exactly zero. 

We shall demonstrate in Sec. IV C that fragmented states allow the occupation of geminals 
that contribute very little to the interaction energy as opposed to condensed states. Thereby, 
the system can lower its energy, once the barrier is high enough. 

B. From condensation to fragmentation 

At barrier heights 13 < A < 24, one finds that the occupation of the second natural orbital 
/N increases continuously from below 1% to almost 50%. The condensate fragments in 
this regime according to the definition of fragmented condensates [23] . In this regime many 
permanents contribute to the wave function and, therefore, we refer to the range of barrier 
heights 13 < A < 24 as the many-body regime. Along with the natural orbital occupations 
the natural geminal occupations change as well. Three natural geminals become occupied 
with increasing barrier height, see Fig. 1. 

In the many-body regime, the upper and lower bounds, Eqs. (5) and (6), on the largest 

(2) (2) 

eigenvalue of the second-order RDM, n\ , no longer restrict n\ to a narrow region. In fact, 

(2) 

n\ takes on a value somewhere in between these bounds. 

This onset of fragmentation manifests itself also in the BMF solution which jumps from 
a GP type permanent \N, 0) in the condensed regime to a fully fragmented solution of the 
form \N/2,N/2). Note that already at barrier heights A > 14 this fragmented solution is 
lower in energy than a GP type permanent. At this barrier height the many-body solution 
is only slightly depleted, see Fig. 1. 

If we compare the natural orbitals in Fig. 2 at barrier height A = 19 with those at A = 13, 
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we note that they look very similar, apart from the fact that the peaks are slightly farther 
apart, and the first natural orbital is closer to zero at the center of the trap. The energies 
of the first two natural orbitals are almost degenerate and the total energy is minimized by 
the occupation of both natural orbitals. Without interactions the system would remain in a 
condensed state, since the single particle energies of the first and the second natural orbitals 
remain separated at any finite barrier height. Note that in the absence of interactions the 
natural orbitals are the eigenfunctions of h. However, as we noted in Sec. IV A, a system of 
repulsively interacting bosons in multi-well traps can lower its energy by occupying several 
natural orbitals, once the barrier is high enough [26-29, 33]. This is precisely the reason for 
the observed onset of fragmentation. 

In the many-body regime the natural geminals are no longer symmetrized products of the 

(2) 

natural orbitals. If we compare a\ (xi,X2) in Fig. 2 at barrier heights A = 13 and A = 19, 
we see that the peaks on the diagonal at Xi = x 2 = ±d/2 decrease, whilst the peaks on the 
off-diagonal at x\ = —x 2 = ±d/2 increase when the barrier is raised. The opposite is true for 
the third natural geminal a^\x±, x 2 ): the off-diagonal maxima at X\ = —x 2 = ±d/2 have 
decreased, whilst the diagonal minima at x\ = x 2 = d/2 are now more negative. On the other 

(2) 

hand, the second natural geminal, a 2 (xi,x 2 ), is still well approximated by a symmetrized 
product of the first and second natural orbital. The behavior of the natural geminals is 
qualitatively different from that displayed by the natural orbitals. In contrast to the natural 
orbitals, the natural geminals do change their shape during the fragmentation transition. 
They only obtain their final forms, when the fragmentation transition is completed, see Fig. 2 
and Sec. IV C. It is possible to understand qualitatively which of the two unoccupied natural 

(2) 

geminals becomes occupied first. In the condensed regime a 2 (xi,x 2 ) is well approximated 
by a state |1, 1) in which one boson occupies the first natural orbital and another one the 
second natural orbital, see Eq. (24). Meanwhile, 0:3 (xi,x 2 ) is well approximated by a state 
1 , 2) , in which two bosons reside in the second natural orbital. Since the single particle 
energies of the first and second natural orbitals are separated, the occupation of the third 
natural geminal costs more energy than that of the second. Hence the second natural geminal 
becomes occupied before the third, see Fig. 1 (bottom). 
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C. Fully fragmented state 



When the central barrier is raised to values A > 24, the two parts of the condensate 
become truly independent. The natural orbital occupations approach = = N/2, 
which reflects the fact that the energies associated with the first and second natural orbitals 
degenerate at infinite barrier heights. The many body wave function can then be adequately 
approximated by a single permanent of the form \N/2, N/2), i.e. with equal numbers of 
particles in the first and the second natural orbitals. Therefore, we refer to barrier heights 
A > 24 as the fully fragmented regime. The additional energy, necessary for the occupation 
of the second natural orbital, is outweighed by a lower interaction energy. Note that this 
final form of the wave function is anticipated by the BMF solution at barrier heights A > 14. 

The natural geminal occupations approach 

nS 2) = N(N/2), nf = = N/2(N/2 - 1) (25) 

in the fully fragmented regime. These are the values that follow from the BMF solution. 

It is only at barrier heights A > 24 that the natural geminals take on their final shapes, 
compare the third and fourth rows of Fig. 2. If we expand the natural geminals in the basis 
of natural orbitals at these barrier heights, we find that 

l«i 2) > = ^= (|2, 0) - |0, 2)) , |«< 2) ) = |1, 1), |c4 2) > = -L (|2, 0) + |0, 2)) (26) 

holds to a very good approximation. The first and third natural geminals have equal con- 
tributions coming from the first and the second natural orbitals. The question arises, why 
their occupations are different, about 50% and 25%, respectively. Subtracting the perma- 
nents |2, 0) and |0, 2) from one another yields a geminal which is localized on the off-diagonal, 
see af\x\, x 2 ) in Fig. 2 at A = 24. Adding the permanents |2, 0) and |0, 2) yields a geminal 
which is localized on the diagonal, see af\x\,x 2 ) in Fig. 2 at A = 24. It is easy to see 
from the shape of the natural geminals in the fourth row of Fig. 2 that the integrals over 
the one-body part in Eq. (8) are approximately the same for each of the natural geminals. 
Given the occupations in Eq. (25), the first natural geminal contributes about one half of 
the one-body energy, whereas the second and the third natural geminal contribute about 
a fourth each. The situation is different for the two-body part of the Hamiltonian. Since 
a 2 2 \x\, x 2 ) and a^\xi,x 2 ) are localized on the diagonal, they do contribute to the interac- 
tion energy. In contrast, a^\xi,x 2 ) is almost zero at coordinate values x\ ~ x 2 and, due 
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to the contact interaction in Eq. (23), it practically does not contribute to the interaction 
energy. At high barriers a fragmented state allows the system to lower its energy through 
the occupation of a natural geminal which is localized on the off-diagonal. 

We would like to make a remark on the validity of the present MCTDHB computation 
for high barriers. For high barriers the whole system can be considered as composed of two 
separate condensates. To describe the depletion of each condensate it would be necessary 
to employ M = 4 orbitals. We use only M = 2 orbitals in the many-body computation and 
cannot describe this depletion. We justify the use of M = 2 orbitals by noting that at A = 
the system is almost fully condensed, and the depletion can be safely neglected, see Fig. I. 
Therefore, we assume that the depletion of each of the two condensates can be neglected, 
when the barrier is very high. This claim is supported by a computation that we carried 
out in the harmonic trap at the same interaction strength Ao = 0.01 for 500 particles. The 
depletion was found to be even less than for N = 1000 particles. 

V. FIRST ORDER CORRELATIONS 
A. General analytical considerations 

We now describe the first-order correlations in an analytical mean-field model for the two 
limiting cases of a condensed and a fully fragmented system. In these cases mean-field theory 
has been shown to be well applicable, see Sec. IV. For our purposes the exact shape of the 
natural orbitals a^\xi) and a^\xi) is unimportant. Consider a normalized one-particle 
function, which is localized at the origin. $(x) may vary in shape, but is always 

assumed to resemble a Gaussian. Similarly, we define translated copies <&i(x) = <&(x + d/2) 
and $20*0 = &(x — d/2) of <&(#), where the previously defined distance d between the minima 
of the potential wells is taken to be large enough to set products of the form ^i(x)^ 2 (x) to 
zero. Since $ is localized in some region around the origin, $1 is localized in a region L to 
the left and $2 in a region R to the right of the origin. 

1. Condensed state 

In the condensed regime, < A < 13, only one natural orbital, a^\xi), is significantly 
occupied. Therefore, we approximate the first-order reduced density operator of the system 

18 



by that of a condensed state \N, 0) 

P$ J0) = N\<*?Wi ) l (27) 

It then follows from Eq. (13) that 

\g\S fi) (^i)\ 2 = l. (28) 

At zero barrier height, the first natural orbital is a Gaussian, broadened by interactions. 
Therefore, we write af\xi) = $(^i), and hence the one-particle density distribution and 
the one-particle momentum distribution are of the form 

A$ >0 >(zil*i) = N \*W\ 2 > ( 29 ) 
A$ >0) (*i|*i) = Nmh)\ 2 . (30) 

Since $(xi) is a broadened Gaussian, its Fourier transform $(fci) is also close to a Gaussian, 
but narrower in comparison to a non-interacting system. The momentum distribution of 
the repulsively interacting system in the harmonic trap is therefore narrower than that of a 
non-interacting system. 

We now turn to the case corresponding to A ^ 13, where the system is still condensed, 
but the first two natural orbitals are spread out over the two wells. We model the natural 
orbitals by 

a?( Xl ) = -±= [^(xO + $ 2 (*i)] , <*?M = ^= ~ $ 2(^i)] • (31) 

In this case one obtains [56]: 

P$o>(sil*i) = yl $ i(^)| 2 + yl^^OI 2 , (32) 
P\N,o)( k ^) = ^V[l + cos(A; 1 d)]|$(A; 1 )| 2 (33) 

for the density and the momentum distribution. We note that the one-particle momen- 
tum distribution displays an oscillatory pattern in momentum space at a period which is 
determined by the separation d of the centers of the two wells. 



2. Fully fragmented state 



In the true many-body regime, 13 < A < 24, where many permanents contribute to the 
wave function, a mean-field model is bound to fail. However, in the fully fragmented regime it 
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is possible to consider the whole system as two separate condensates, and hence a mean-field 
description is again applicable. Therefore, we now turn to the case corresponding to A > 24, 
where the system is fully fragmented and the many-body state is given by \N/2, N/2). The 
first-order density operator then reads: 

p\$ /2 , Nm = + yla^X^I- (34) 

Since the natural orbitals remain qualitatively unchanged during the fragmentation transi- 
tion, we approximate a^\xi) and a^\x\) by Eqs. (31) and obtain for the density and the 
momentum distribution [56]: 

P%2,n/2)(^W) = ylMzi)| 2 + yl$2(xi)| 2 , (35) 
$/2,W*il*i) = N \*fa)\ 2 . (3 6 ) 

We note that the one-particle momentum distribution of independent condensates does not 
contain an oscillatory component and is identical to the momentum distribution of a single, 
localized condensate of N particles within this model, see Eq. (30). For the normalized 
first-order correlation function one finds 



1 if xi, x[ G L or x\, x[ G R, 
otherwise. 



I (1) // n,2 ) "-1,-1^-^-1,-1^^, 

\9\N/2,N/2)y X l' ~~ \ \ 61 ) 



Whereas the state \N, 0) is fully first-order coherent, the fragmented state \N/2, N/2) is 
only first-order coherent in a restricted and generally disconnected region. Each of the 
two condensates is first-order coherent, but the mutual coherence which is present in the 
condensed regime is lost. 



B. Numerical results 



We now turn to the discussion of first-order correlations. In particular, we are interested in 
effects that are due to the true many-body nature of the wave function. Along with our many- 
body results we plot the corresponding results of the BMF solution. From the discussion 
in Sec. IV it is clear that we expect many-body effects to occur during the fragmentation 
transition at barrier heights 13 < A < 24. In the condensed and in the fully fragmented 
regime we expect that the many-body results are well approximated by those of the BMF 
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solution. In these cases we can understand the structure of the correlation functions on the 
basis of the analytical mean-field model of Sec. V A. 

The first column of Fig. 3 shows the one-particle density distribution p^(xi\xi) of the 
many-body solution (blue line) and that of the BMF solution (red line with triangles) at the 
barrier heights A = 0, 13, 19, 24, from top to bottom. It is remarkable that the one-particle 
densities obtained from either the many-body wave function or the BMF solution give results 
that cannot be distinguished from one another at any barrier height. 

In a purely harmonic trap, A = 0, the one-particle density takes on the form of an 
interaction-broadened Gaussian. At higher barriers, the density splits into two parts that 
are localized in each of the wells. At A = 13, the one-particle density has developed two 
separated peaks. Note that the system is still in the condensed regime at this barrier height 
and must be considered a single condensate, despite the spatial separation between the two 
peaks. 

When the central barrier is raised further to values 13 < A < 24 the condensate fragments, 
see Fig. 1. At a barrier height of A = 19 the system is halfway on its way from a condensed 
to a fully fragmented condensate. Many permanents contribute to the many-body wave 
function and one may wonder how this transition manifests itself in observable quantities. 
However, apart from a small shift of the center of the two peaks and a reduction at the 
center of the trap, p^(xi\xi) remains largely unaffected by this transition. If the barrier 
is raised further to A = 24, the fragmentation transition is largely completed. Also during 
the transition from a true many-body state to a fully fragmented state there is no visible 
indication of this transition in the one-particle density. 

The second column of Fig. 3 shows the one-particle momentum distribution p^\ki\ki) 
at the same barrier heights as before. At A = 0, the one-particle momentum distribution 
is given by a squeezed Gaussian, in agreement with Eq. (30). At A = 13 the one-particle 
momentum distribution has developed an oscillatory pattern, typical of a single condensate 
spread out over two wells. The structure of p^(ki\ki) is well reproduced by Eq. (33) of the 
analytical mean-field model. Up to this barrier height the BMF solution is almost identical 
to the many-body wave function, and therefore the respective momentum distributions are 
indistinguishable, see the two upper panels in the second column of Fig. 3. 

When the system enters the many-body regime, 13 < A < 24, the momentum distribution 
of the many-body solution deforms to a Gaussian-like envelope, modulated by an oscillatory 
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part. The BMF momentum distribution, on the other hand, already takes on the form 
characteristic of two separate condensates. It agrees with the prediction of Eq. (36), which 
is clearly different from the many-body result. This merely reflects the fact that the exact 
wave function is inaccessible to mean-field methods in the many-body regime. 

When the state becomes fully fragmented at A = 24, the many-body momentum dis- 
tribution and the BMF momentum distribution become indistinguishable again, consistent 
with an explanation in terms of two independent condensates, see Eq. (36). Compared to 
p^(ki\ki) at A = 0, the momentum distribution is broader at A = 24, because the density 
distribution in each of the two wells is narrower than that in the harmonic trap. 

The third column of Fig. 3 shows the absolute value squared of the normalized first-order 
correlation function \g^(x[, Xi)\ 2 of the many-body solution only. Here and in all following 
graphs of correlation functions we restrict the plotted region by a simple rule. To avoid 
analyzing correlations in regions of space where the density is essentially zero, we plot the 
respective correlation function only in regions where the density is larger than 1% of the 
maximum value of the density in the entire space. We apply the same rule also in momentum 
space. 

At zero barrier height \g^(x[, xi)\ 2 is very close to one in the region where the density 
is localized. The system is first-order coherent to a very good approximation and the mean- 
field formula Eq. (28) applies. As the barrier is raised to A = 13 the coherence between the 
two peaks, e.g. at x 1 = —x[, is slightly decreased, while the coherence within each of the 
peaks is preserved. Note that the density at the center of the trap is already below 1% of the 
maximal value in this case. Despite this separation the system remains largely condensed, 
but deviations from Eq. (28) are visible. If the barrier is raised further to A = 19, the 
coherence of the system on the off-diagonal decreases quickly. Although the bosons in each 
well remain coherent among each other, the overall system is only partially coherent. At 
barrier heights A > 24, the coherence between the two wells is entirely lost. This is also the 
scenario that the BMF solution anticipates, see Eq. (37). 

It is remarkable that not only the density, but also the momentum distribution obtained 
within mean-field theory agree so well with the many-body result, when the system is not in 
the true many-body regime. This would not be the case if we had restricted the mean-field 
approach to the GP equation, as we shall show now. 

Up to barrier heights A = 13 the many-body system is condensed, and the BMF solution 
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coincides with the GP solution. The BMF, and therefore also the GP solution provide good 
approximations to the interacting many-body system. Above A = 13 the results obtained 
with the GP mean-field become qualitatively wrong as the barrier is raised. To illustrate 
this point, we plot the GP results corresponding to those of Fig. 3 at barrier heights A = 19 
(top) and A = 24 (bottom) in Fig. 4. A comparison of the respective one-particle densities, 
shown in the first column of Fig. 3 and Fig. 4, reveals no visible difference. The GP mean- 
field reproduces the density distribution at all barrier heights correctly. However, the GP 
solution fails at the description of the momentum distribution and the normalized first order 
correlation function, compare the second and third columns of Figs. 3 and 4 at the same 
barrier heights. The reason for the failure of the GP mean-field is the assumption that 
all bosons occupy the same orbital. It is by construction incapable to describe fragmented 
condensates. 

VI. SECOND ORDER CORRELATIONS 
A. General analytical considerations 

In this subsection we extend the analytical mean-field model of Sec. VA to describe 
second-order correlations. 

1. Condensed state 

We found in Sec. IV A that only one natural geminal is significantly occupied in the 
condensed regime, where the many-body state is approximately given by a single permanent 
in which all bosons occupy the same single particle state. Therefore, we approximate the 
second order reduced density operator in the condensed regime by that of the state | N, 0): 



natural orbital a\ ' . For the condensed state \N, 0) one finds that up to corrections of order 




(38) 



where a\ (x±, x?) 



a^\xi)cx[\x2) is the permanent in which two bosons reside in the first 
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0(1/N) the state is second-order coherent: 

g$ 0) (x 1 ,x 2 ,x 1 ,x 2 ) = 1 - — , (39) 

gfl^hMMM) = i-^- (40) 

Thus, there are practically no two-body correlations if AT ^> 1. At zero barrier height the 
first natural orbital takes on the shape of a broadened Gaussian, a^\xi) = where 
$(x) is defined in Sec. VA. The first natural geminal then reads oc[\xi,x 2 ) = $(xi)$(x 2 ). 
It follows that the two-particle density and the two-particle momentum distribution factorize 
up to corrections of order 0(1/ N) into products of the respective one-particle distributions: 

p\^ 0) (x u x 2 \x u x 2 ) = N(N-l)m Xl )\ 2 mx 2 )\ 2 , (41) 
pfl G) {k u k 2 \kiM) = iV(iV -1)|$(A; 1 )| 2 |$(A; 2 )| 2 . (42) 

At the barrier height A = 13, the system is condensed but spread out over the two wells. 
Then, using Eqs. (31) to approximate we find 

pf^ 0) (xi,x 2 |xi,x 2 ) = Ni < N ~ l A [|$ 1 ( 3 ; 1 )'l» 1 (3; 2 )| 2 + \$ 1 (x 1 )$ 2 (x 2 )\ 2 + 

|$ 2 (zi)* i{x 2 )\ 2 + |^ 2 (^i)^2(^ 2 )| 2 ] , (43) 
pfl^ihMkiM) = N{N -l)[l + cos(M)][l + cos(M)]|<f(A;i)<f(A; 2 )| 2 (44) 

for the two-particle density and the two-particle momentum distribution. Apart from a 
correction of order 0(1/N), the two-particle density and the two-particle momentum distri- 
bution are again products of the respective one-particle distributions. 

2. Fully fragmented state 

In Sec. IV C we found that three natural geminals are occupied in the fully fragmented 
regime, see Eq. (25). The occupations of Eq. (25) hold exactly for a state of the form 
\N/2,N/2). Therefore, we approximate the second-order reduced density operator in the 
fully fragmented regime by that of the state \N/2, N/2): 

= N J l«i 2) ><«S 2, l + T (t " }a?) ><Q ' f 1 + 7 (f - |af)<Qf 1 (45) 

where the natural geminals \af^) are given by Eq. (26). In contrast to the condensed state, 
the normalized second-order correlation function of the fully fragmented state has a more 
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complicated structure due to the different terms contributing to Eq. (45). We approximate 
the natural geminals using Eqs. (31) and find 

~ 0<J>i(*i)<J> 2 (* 2 )| 2 + Wsi)* i(^)| 2 ] , (46) 
P% 2 , N/2) (h,k 2 \ kl ,k 2 ) = N(N-l) (l + jV ^ 1 COS[(A:i 2 " A:2M] ) l*(*i)*(k)l 2 ( 4 7) 

for the two-particle density and the two-particle momentum distribution. This represen- 
tation allows us to identify the first two terms in Eq. (46) as contributions coming from 
two separate condensates of N/2 bosons each, with condensate wave functions 3>i(xi) and 
$2(^1) ■ The third term in Eq. (46) is due to the fact that the bosons in the two separated 
condensates are identical particles. For the normalized second-order correlation function one 
finds: 

(2) I 1 - if Xi, x 2 G L or xi, x 2 G R 

9\N/2,N/2)( X 1> X *> X 1> X 2) = ) > ( 48 ) 

I 1 otherwise 

which mimics a high degree of second-order coherence. However, when is evaluated on 
the diagonal in momentum space, one finds 

$W*^.*»> = (l 4) (l + J^H^l^M) , (49) 

which displays an oscillatory behavior and deviates significantly from a uniform value of 
one. Hence the system is clearly not coherent, see Sec. II. Therefore, a description of second 
order correlations in terms of p^ 2 \xi,x 2 \xi,x 2 ) and g^ 2 \xi,x 2 ,xi,x 2 ) is incomplete, and 
p^(k 1: k 2 \k 1: k 2 ) and g^ 2 \k-i,k 2 ,ki,k 2 ) have to be taken into account. Although this may 
seem obvious in the present case of a fully fragmented state, this reduction of coherence is 
more intricate in a state which is only partially fragmented, see following subsection. 



B. Numerical results 



In this subsection we discuss the second-order correlations of the many-body solution. 
We compare the results to those of the BMF solution. When mean-field theory gives a good 
approximation to the many-body results, we also compare with the analytical mean-field 
model of Sec. VI A. 
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The first two columns of Fig. 5 show the two-particle density p^(x 1: x 2 \x 1: x 2 ) of the 
many-body (left) and BMF (right) solutions at barrier heights A = 0, 13, 19, 24, from top to 
bottom. At zero barrier height p^ 2 \xi,x 2 \xi,x 2 ) is localized at the center of the trap. The 
two-particle density factorizes approximately into a product of the one-particle densities: 
p( 2 \xi, x 2 \xi, x 2 ) ~ p^(xi\xi)p^(x 2 \x 2 ). This remains true up to barrier heights A = 13, 
where the condensate is spread out over the two wells. The BMF result approximates the 
many-body result well in the condensed regime, and the structure of p^ 2 \x\, x 2 \xi, x 2 ) is 
that of Eqs. (41) and (43) at barrier heights A = and A = 13, respectively. 

When the barrier is raised further to A = 19, the system fragments. Many permanents 
contribute to the wave function in this regime and there is no simple formula that relates 
the occupations of the natural orbitals to the two-particle density. Similar to the one- 
particle density, described in Sec. VB, the two-particle density seems to take no notice of 
the transition from a single to a fragmented condensate. It remains practically unchanged 
during the transition, apart from a slight shift of the peaks away from each other as the 
barrier is raised. 

At even higher barriers, A > 24, the many-body state becomes fully fragmented and 
the wave function approaches \N/2, N/2). In this limit it is again possible to describe the 
two-particle density on a mean-field level. Therefore, the analytical results of Sec. VI A for 
the fully fragmented state should apply. In fact, the structure of p^(xi,x 2 \xi,x 2 ) in the 
fully fragmented regime is that predicted by Eq. (46). 

The two-particle density of the condensed state just below the fragmentation transition 
and of the fully fragmented state above the fragmentation transition cannot be distinguished. 
It is easily verified, that Eqs. (43) and (46) give rise to the same two-particle density profile 
up to corrections of order 0{1/N). 

In contrast, the fragmentation transition is clearly visible in the two-particle momen- 
tum distribution. In the third and fourth columns of Fig. 5 the two-particle momentum 
distribution p^ 2 \ki, k 2 \k\, k 2 ) of the many-body (left) and BMF (right) wave function are 
shown. 

In the condensed regime the two-particle momentum distribution is approximately given 
by the product of one-particle momentum distributions of a single condensate. This agrees 
with the analytical predictions of Eq. (42) at barrier height A = and Eq. (44) at A = 13. 
The mean-field picture is appropriate here. 
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In the many-body regime the two-particle momentum distribution p( 2 \ki, k 2 \k 1 , k 2 ) con- 
tains contributions from many permanents. The resulting p( 2 \k\, k 2 \ki, k 2 ) has a structure 
that lies somewhat in between the two results, Eqs. (44) and (47), obtained within the 
analytical mean-field model. The BMF solution is fully fragmented and does not provide 
an accurate approximation to the many-body two-particle momentum distribution in this 
regime, see Fig. 5, third and fourth columns in the third row from above. 

When the barrier is raised to A = 24, the many-body state becomes fully fragmented and 
the mean-field picture is again applicable. The pattern of a single coherent condensate has 
now vanished completely in favor of a pattern characteristic of two separate condensates. 
The pattern agrees well with the structure predicted by Eq. (47). 

Similar to our results on first order correlations, discussed in Sec. VB, the fragmentation 
transition shows up in the two-particle momentum distribution, but not in the two-particle 
density. While this behavior is predictable in the limiting cases of a condensed and a fully 
fragmented state, it is necessary to solve the many-body problem to determine the limits 
of such mean-field approximations. Particularly the behavior in between the two mean-field 
limits is only accessible to many-body approaches. 

We will now address the second-order coherence of the system. The first two 
columns of Fig. 6 show the diagonal of the normalized second order correlation function 
g( 2 \xi,x 2 ,Xi,x 2 ) of the many-body (left) and the BMF (right) solutions. Note the scale! 
The Eqs. (39) and (48) of the analytical mean-field model of Sec. VIA predict very small 
correlations in the two-particle density of the condensed and the fragmented state. This is 
confirmed in the first column of Fig. 6. In the condensed regime at zero barrier height the 
effects of the depletion of the condensate on g( 2 \xi,x 2 ,Xi,x 2 ) are visible. Almost no two- 
particle density correlations are present. This is equally true in the case of a single condensate 
spread out over the two wells and also in the many-body regime. Above the fragmentation 
transition, the present computation of the many-body solution cannot describe effects on 
g^> that are due to the depletion of the condensate. However, since depletion effects are 
negligible in the harmonic trap, we are reassured that they are also negligible in the fully 
fragmented regime, see Sec. IV C. The BMF solution predicts almost identical two-body 
density correlations, see second column of Fig. 6. 

On the basis of g( 2 \xi, x 2 , x±, x 2 ) alone, the many-body state appears to be second-order 
coherent at all barrier heights. A high degree of second-order coherence requires Eq. (14) 
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to hold to a very good approximation for p = 1 and p — 2. This in turn requires the 
largest eigenvalues of the first- and second-order RDM to be m N and N(N — 1), 

respectively. We have already demonstrated in Sec. IV that these conditions are only satisfied 
in the condensed regime. Therefore, it is obviously tempting, but wrong to conclude from 
g^(xi,X2,xi,X2) ~ 1 that the system is second-order coherent. This misconception is due 
to the fact that g^(xi,X2,Xi,X2) only samples a small part of the first and second-order 
RDMs of the system. 

So, how does the decrease of coherence manifest itself in second order correlation 
functions? For second-order coherence to be present, at least approximately, also 
g( 2 \k\, &2, &4, k 2 ) has to be close to one. The third and fourth column of Fig. 6 show 
g^(ki,k2,ki,k 2 ) of the many-body (left) and BMF (right) solution. At zero barrier height 
the system is indeed highly second-order coherent since only one natural orbital is signifi- 
cantly occupied. Not only g( 2 \xi, x 2 , Xi, x 2 ), but also g( 2 \k\, k 2 , ki, k 2 ) is very close to one 
here. However, at A = 13 when the many-body state is still condensed, g( 2 \ki, k 2 , fci, k 2 ) 
starts to develop a structure. 

When the barrier is raised to values above A = 13, the structure becomes more and 
more pronounced. In the many-body regime at A = 19, we find that g( 2 \k±, k 2 , k±, k 2 ) has 
a complicated behavior and deviates significantly from values close to one, thereby proving 
that strong correlations are present. Note that the interaction between the particles is weak 
and that the strong correlations are due to the transition from a single to a fragmented 
condensate. This transition is in turn induced by a change of the shape of external potential. 
Varying the shape of the external potential therefore provides a means to introduce strong 
correlations between the particles. The strongest correlations (black spots in the third panel 
of the third row of Fig. 6) occur at those values where the two-body momentum distribution 
has local minima. At the values of k\ and k 2 , where the strongest correlations occur, the 
one-body and the two-body momentum distributions are clearly distinct from zero, see 
third panel in the middle column of Fig. 3 and the third panel in the third row of Fig. 5. 
Experiments that measure g^^ki, k 2 , ki, k 2 ) in ultracold quantum gases have been carried 
out recently, see e.g. [57]. An experiment that measures g( 2 \k\, k 2 , k±, k 2 ) would find strong 
two-particle momentum correlations at high barriers. 

When the system becomes fully fragmented at barrier heights A > 24 the structure of 
g( 2 \ki, k 2 , ki, k 2 ) becomes more regular again. The amplitude of the correlations is smaller 
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than in the many-body regime, and the correlations between different momenta are mod- 
ulated by a single oscillatory structure. This structure can be well understood within the 
analytical mean-field model of Sec. VIA. The oscillatory part of g^ 2 \ki, k 2 , fci, k 2 ) is deter- 
mined by the difference of the wave vectors multiplied by the distance between the wells, 
see Eq. (49). 

Hence, we find that only in the condensed regime the system is second-order coherent 
despite the fact that g^(x 1: x 2 , x 1: x 2 ) ~ 1 at all barrier heights. This merely reflects the 
fact that g ( - 2 \xi,x 2 ,Xi,x 2 ) is only the diagonal of g^(x[, x' 2 , Xi, x 2 ). On the other hand, 
g( 2 \ki, k 2 , ki, k 2 ) depends on all values of of p^(xi, x 2 \x\, x 2 ) and provides complementary 
information about the coherence of the state. A description of second-order coherence in 
terms of g^ 2 \x\, x 2 , x±, x 2 ) alone is therefore incomplete. 

The corresponding results of the BMF solution agree well with those of the many-body 
solution as long as the system is not in the many-body regime at intermediate barrier heights. 
In the many-body regime the BMF result is inaccurate, but it anticipates the final form of 
g( 2 \xi, x 2 , xi, x 2 ) in the fragmented regime. 

VII. CONCLUSIONS 

In this work we have investigated first- and second-order correlations of trapped inter- 
acting bosons. For illustration purposes we have investigated the ground state of N = 1000 
weakly interacting bosons in a one-dimensional double-well trap geometry at various bar- 
rier heights on a many-body level. We have obtained the many-body results by solving 
the many-body Schrodinger equation with the recently developed MCTDHB method. This 
allowed us to compute from first principles the natural orbitals and the natural geminals 
of a large interacting many-body system, together with their occupation numbers. To our 
knowledge this is the first computation of the natural geminals of an interacting many-body 
system of this size. 

Depending on the height of the double-well barrier we found that there are three different 
parameter regimes. At low barriers the ground state is condensed and the many-body wave 
function is well approximated by a single permanent of the form \N, 0). At high barriers the 
ground state becomes fully fragmented and can be well approximated by a single permanent 
of the form \N/2, N/2). At intermediate barrier heights, where the transition from a single to 
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a fragmented condensate occurs, the ground state becomes a true many-body wave function 
to which many permanents contribute. We have demonstrated that the transition to a 
fragmented state results in the occupation of a natural geminal that contributes very little 
to the interaction energy. The overall energy of the system can be lowered by the occupation 
of such a geminal, and the ground state becomes fragmented. 

We have shown how the transition from a condensed to a fully fragmented ground state 
manifests itself in the one- and two-particle momentum distributions. However, the tran- 
sition is not captured by the one- and two-particle density distributions, not even in the 
many-body regime. 

In order to determine the coherence of the state during the fragmentation transition, 
we have computed the first- and second-order normalized correlation functions g( 1 \x' 1 ,xi), 
g( 2 \x\, X2, xi,x 2 ) and g( 2 \ki, k 2 , fci, k 2 ). In the condensed regime, a high degree of coherence 
is indeed present in the ground state wave function. First and second order correlations 
were found to be negligible at the interaction strength and particle number chosen for our 
computation. However, with increasing barrier height correlations between the momenta of 
the particles build up. These correlations were found to be very strong in the many-body 
regime at intermediate barrier heights. The ground state at high barriers was found to be 
correlated, but not as strongly as the ground state at intermediate barrier heights. 

While the transition from a virtually uncorrelated state to a correlated one is clearly visi- 
ble in g^(x[, Xi) and g( 2 \ki, k 2 , k±, k 2 ), the transition hardly shows up in g( 2 \xi, x 2 , Xi, x 2 ). 
A description of second-order coherence in terms of g^(xi,x 2 ,Xi,x 2 ) alone is, therefore, 
incomplete and can lead to wrong predictions. 

For comparison we have computed results based on (i) the best approximation of the 
many-body wave function within mean-field theory, the BMF wave function, and (ii) the 
Gross-Pitaevskii solution. We found that the GP wave function is identical to the BMF 
solution up to some barrier height. However, once the true many-body solution starts to 
fragment the BMF wave function is no longer given by a GP type permanent \N, 0), but 
rather by a fragmented state of the form | AT/2, N/2). In the true many-body regime neither 
the GP, nor the BMF solution provide an adequate approximation to the many-body wave 
function, and the predicted correlations are inaccurate. 

While the GP mean-field is only accurate at low barrier heights, the BMF solution pro- 
vides a very good approximation to the true many-body wave function at low and high 
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barriers. We have shown that the GP mean-field predicts qualitatively wrong results at high 
barriers. The BMF only fails at intermediate barrier heights where the true many-body 
wave function becomes a superposition of many permanents. Such many-body effects can, 
by construction, not be captured by mean-field methods. 

In the mean-field regimes at high and low barriers we have provided an analytical mean- 
field model that allows us to understand the general structure of the computed correlation 
functions. 

Our work sheds new light on the first- and second-order correlation functions of interacting 
many-body systems. The variation of the shape of the trapping potential allows one to 
change the physics of the system from mean-field to strongly correlated many-body physics. 
Particularly, the many-body regime in between the condensed and the fully fragmented 
regimes has shown to be very rich and promises exciting results for experiments to come. 
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APPENDIX A: p-PARTICLE MOMENTUM DISTRIBUTION 

The relation of the p-particle momentum distribution to the p-th order RDM is shown. 
The p-particle RDM is related to the p-particle momentum distribution by 

p^(k 1 , . . . , k^lkx, . . . , k p ; t) = f dPvdPv' k ^- r ;) p (p)( ri) . . . , r p |r;, . . . , r; ; t). 

(Al) 

The change of variables Rj = ^rp-, Sj = r; — for i — 1, . . . ,p in Eq. (Al) leads to 

p^(k 1 ,...,k p |k 1 ,...,k p ;t) = yVse-^ k < s <7 {p) (si,...,s P ;t), (A2) 

where 

7 W(s 1 ,...,s p ;f) = J ^Rp(rt(R 1 + |,...,R p + || Rl -|,...,R p -| ; t) (A3) 

is the average of the value of p^(ri, . . . , r p \r' 1: . . . , r' p ; t) at distances between and r-. 
From Eqs. (A2-A3) it is clear that the p-particle momentum distribution at large momenta 
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is determined by the behavior of p^ p \r 1: . . . ,r p \r' 1: . . . ,r' p ;t) at short distances, whereas at 
low momenta the off-diagonal long range behavior of p^(ri, . . . ,r p \r[, . . . ,r' p ;t) contributes 
the major part. 
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FIG. 1: (color online). Energy per particle, natural orbital and natural geminal occupations of 
the ground state of N = 1000 bosons at Ao = 0.01 in a harmonic trap with a central barrier. 
Shown is the dependence on the barrier height. Top: energy per particle E/N of the many-body 
solution. Inset: energy difference per particle between the best mean- field and the GP solution, 
(E 'bm f — Egp) l 'N (triangles), and between the many-body and GP solution, (Emctdhb—Egp)/N 
(circles). Middle: the eigenvalues and of the first-order RDM p^ l \x\\xi). The ground 

(2) (2) (2) 

state fragments with increasing barrier height. Bottom: the eigenvalues n\ ,nr> and of the 
second-order RDM p( 2 >(xi,X2\x' 1 ,x' 2 ). The dashed lines in the middle and bottom panel indicate 
upper and lower bounds on the largest eigenvalue of the first- and second-order RDMs. See text 
for details. The quantities shown are dimensionless. 
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FIG. 2: (color online). Natural orbitals and geminals at different barrier heights. First column: the 
natural orbitals a^\xi) (red line) and a^\x±) (blue line) of the many-body solution at different 
barrier heights ^4 = 0, 13, 19, 24, from top to bottom. The trapping potential is shown as a black 
line in the first column. The state of the system changes from condensed to fragmented between 

(2) (2) 

A = 13 and A = 24. Second to fourth columns: natural geminals a\ (xi,X2), «2 (xi,%2) an d 

(2) 

ag (xi,X2) from left to right at the same barrier heights as above. While the natural orbitals 
remain qualitatively unchanged during the fragmentation transition, the natural geminals take 
on their final shapes only when the system becomes fully fragmented. At low barrier heights 
only one natural geminal is occupied. At high barriers three natural geminals are occupied, see 
Fig. 1. The total energy is minimized by the occupation of a natural geminal that contributes 
practically nothing to the interaction energy. See text for more details. The quantities shown are 
dimensionless. 
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FIG. 3: (color online). Density distribution, momentum distribution and first-order coherence. 
The first two columns show the one-particle density p^\xi\xi)/N and the one-particle momentum 
distribution p^\ki\ki)/N of the many-body solution (blue line) and of the BMF solution (red 
line with triangles), respectively. From top to bottom the height of the central barrier is A = 
0, 13, 19, 24. The BMF result agrees well with the many-body result for a large range of barrier 
heights. At A = 19, in the many-body regime, deviations are visible in the momentum distribution. 
See text for details. The third column shows the absolute value squared of the normalized first-order 
correlation function Ig^'^x^, x\)\ 2 at the same barrier heights. An initially coherent condensate 
splits into two separate condensates which are no longer mutually coherent. Only the coherence 
within each of the wells is preserved. The quantities shown are dimensionless. 
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FIG. 4: (color online). Density distribution, momentum distribution and first-order coherence 
obtained by using the GP equation for high barriers. The first two columns show the GP one- 
particle density p^ 1 \x\\xi) /N (left) and the GP one-particle momentum distribution p^ l \k\\ki)/N 
(middle) at barrier heights A = 19 and A = 24 (solid green lines). In the first column the trapping 
potential is also shown (solid black line). The GP equation models the density well, but fails at 
the computation of the momentum distribution, compare with Fig. 3. The third column shows the 
absolute value squared of the normalized first-order correlation function ^(^(x^, xi)| 2 computed 
with the GP equation at the same barrier heights. The normalized first-order correlation function is 
incorrectly described by the solution of the GP equation. The quantities shown are dimensionless. 
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FIG. 5: (color online). Two-particle density and two-particle momentum distribution at differ- 
ent barrier heights. The first two columns (from left to right) show the two-particle density 
p( 2 \x\, X2\xi, X2)/N(N — 1) of the many-body (left) and BMF (right) wave function for the barrier 
heights A = 0, 13, 19, 24, from top to bottom. At low barrier heights (A = 0, 13) the system is con- 
densed, and the two-particle RDM factorizes into a product of the one-particle RDMs. At higher 
barriers (A = 19,24), the system fragments and the two-particle RDM does not factorize into a 
product of the one-particle RDMs. The two-particle density at high barriers looks very similar to 
that of the condensed state at A = 13. The fragmentation transition is not visible in the two- 
particle density. The results of the many-body and BMF wave function cannot be distinguished 
at any barrier height. The third and fourth column show the two-particle momentum distribution 
p( 2 \k\, k2\k\, k2)/N(N — 1) of the many-body (left) and BMF (right) solution at the same barrier 
heights as above. The transition from a condensed state to a fragmented state is clearly visible. 
At A = 19 the BMF solution does not reproduce the many-body results. The system is in a 
true many-body state, inaccessible to mean-field methods. At even higher barriers A > 24 the 
system fully fragments, and a mean-field description is applicable again. The quantities shown are 
dimensionless. 
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FIG. 6: (color online). Second order coherence at different barrier heights. The first two columns 
(from left to right) show the diagonal of the normalized second-order correlation function in real 
space g( 2 \x\, X2, xi, x%) of the many-body (left) and BMF (right) solution at barrier heights A = 
0, 13, 19, 24, from top to bottom. g^ 2 \x\, X2, x±, X2) is very close to one at all barrier heights. Note 
the scale! The system seems to be second-order coherent and the results of the many-body and 
BMF solution agree well with each other. The third and fourth column depict the diagonal of the 
normalized second-order correlation function in momentum space g@)(ki, &2, k\, A^) of the many- 
body (left) and BMF (right) solution at the same barrier heights. The fragmentation transition is 
clearly visible between A = 13 and A = 24. At A = 19 there are strong many-body correlations 
between the momenta (local maxima in black color) and g <y2 \k\, &2, ki, k%) exhibits a complicated 
pattern, see text for more details. A mean-field description fails here. In the limit of high barriers, 
A > 24, the correlations of the many-body state become again describable by those of the BMF 
solution. The quantities shown are dimensionless. 
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